Nonlinear excitations in arrays of Bose-Einstein condensates 
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The dynamics of localized excitations in array of Bose-Einstein condensates is investigated in the 
framework of the nonlinear lattice theory. The existence of temporarily stable ground states dis- 
playing an atomic population distributions localized on very few lattice sites (intrinsic localized 
modes), as well as, of atomic population distributions involving many lattice sites (envelope soli- 
tons), is studied both numerically and analytically. The origin and properties of these modes are 
shown to be inherently connected with the interplay between macroscopic quantum tunnelling and 
nonlinearity induced self-trapping of atoms in coupled BECs. The phenomenon of Bloch oscillations 
of these excitations is studied both for zero and non zero backgrounds. We find that in a definite 
range of parameters, homogeneous distributions can become modulationally unstable. We also show 
that bright solitons and excitations of shock wave type can exist in BEC arrays even in the case 
of positive scattering length. Finally, we argue that BEC array with negative scattering length in 
presence of linear potentials can display collapse. 
PACS numbers: 42.65.-k, 42.50. Ar,42.81.Dp 



I. INTRODUCTION 

The realization of Bose-Einstein condensates (BEC) 
and the observation of quantum interference phenomena 
between two coupled condensates has opened a new fas- 
cinating field in physics, with the perspective of getting a 
better understanding of the complicate behavior of quan- 
tum many body systems and with the hope of realizing 
novel concrete applications of quantum mechanics such 
as atom interferometers and atom lasers. The existence 
of periodic localized oscillations in the relative atomic 
population of one of two coupled condensates (Bloch os- 
cillations), was first predicted in Refs. The close 
analogy existing with the Josephson effect was also em- 
phasized, the common origin of these phenomena being 
the temporal interference of two macroscopic quantum 
states, leading to current oscillations in Josephson junc- 
tions and to oscillations of the atomic population in 
coupled BECs @. 

Although the generalization of these results to the case 
of three coupled BECs was considered , few theoretical 
investigations exist on interference phenomena in arrays 
of coupled BECs 0. In a recent paper the prob- 
lem of Bloch oscillations of bright solitons was investi- 
gated in terms of a tight-binding model for BEC arrays 
with positive scattering length. It is a fortunate situation 
that the equations arising in this case, formally coincide 
with those studied in the theory of nonlinear lattices || . 
One can indeed transpose the knowledge gained in these 
fields, to the field of BEC arrays. Thus, for example, 
one can expect that the types of elementary excitations 
which can arise in BEC arrays, as well as their dynam- 
ical behavior, strongly depend on the ratio between the 



coupling of adjacent condensates and the nonlinearity in- 
duced by the inter atomic interactions. Keeping the same 
terminology used in nonlinear lattices theory, one can say 
that for weak coupling and strong nonlinearity, intrinsic 
localized modes (ILM), i.e. matter excitations localized 
on few lattice sites, should exist. These modes could 
arise in experiments like the ones reported in Ref. [Tc| ] 
in which millions of atoms were trapped in an almost 
one-dimensional ID geometry ("cigar" shape geometry), 
with a small tunnelling term between condensates. 

On the other hand, when the coupling constant is com- 
parable with the self-trapping interactions, small ampli- 
tude excitations of large size (in comparison with the 
lattice constant), should arise. These excitations can 
be identified as lattice envelope solitons (ES) and could 
be observable in macroscopic quantum interference ex- 
periments like the ones reported in Ref. pd[ | in which, 
vertical arrays of "pancake" -like BECs, each containing 
thousands of 87 Rb atoms, were coupled through the grav- 
itational field. In this case the interference phenomenon 
shows up as Bloch oscillations of the ES describing the 
atomic population along the array (ES dynamical local- 
ization). Bloch oscillations of ES type are also possi- 
ble in horizontal optical lattices induced by two counter- 
propagating laser beams with a frequency detuning vary- 
ing linearly in time (this produces a linear potential 
on the BEC similar to the gravitational field of the ver- 
tical traps). In the following we shall call both types 
of excitations (i.e. ILM or ES) as discrete matter soli- 
tons (DMS) whenever it will be clear from the context to 
which type we refer [T^ ] . 

The aim of this paper is twofold. From one side, we 
wish to make a bridge from the theory of nonlinear lat- 
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tices to the field of BEC arrays. This will allow to get 
a better understanding of the types of excitations which 
can arise in ID coupled BECs as the physical parame- 
ters are varied. From the other side, we wish to expand 
both theoretically and numerically, on the phenomenon 
of Bloch oscillations in the case of ID BEC arrays in 
presence of uniform backgrounds (with the term "back- 
ground" here we mean the presence of a uniform distribu- 
tion of atomic population along the array). In particular, 
we discuss: 

i) the existence of ILM in ID BEC arrays with posi- 
tive scattering length; 

ii) the existence of ES in BEC arrays with positive 
scattering lengths; 

iii) Bloch oscillations of soliton-like excitations on top 
of uniform backgrounds in periodic BECs with pos- 
itive scattering lengths; 

iv) the dynamics of BEC arrays with negative scatter- 
ing length on zero backgrounds. 

As to point i) we remark that in contrast with the 
solitons and breathers discussed of Ref . 1 1] which are ex- 
tended on many lattice sites jL4|, the bright ILM dis- 
cussed here are highly discrete, i.e. localized on few lat- 
tice sites. Since ILM have not been observed in real BEC 
arrays, this is an implicit call for experimental work to be 
done in this direction (for experimental evidence of ILM 
existence in other physical contexts see Point ii) 

generalizes the results of || to the case of non-zero back- 
grounds, while in iii) we show that modulational insta- 
bility can arise as the slope of the linear potential (grav- 
itational field for vertical arrays) is increased. Moreover, 
we find that for certain parameter values, shock wave 
excitations can develop. In point vi) we discuss Bloch 
oscillations of bright solitons both analytically and nu- 
merically. All these results will have a close analogy with 
those derived in the theory of nonlinear lattices on in- 
trinsic localized modes p^l9| , envelope solito ns [p(i| (22[ . 
shock waves 21,23}], and Bloch oscillations |2^ , |25|] . 

In Section II the 
In Section III 



II. THE MODEL 



The paper is organized as follows, 
model of a ID array of BEC is derived, 
we discuss the existence of intrinsic localized BEC soli- 
tons, both of even and odd symmetry. In Section IV we 
use a multi-scale expansion to discuss the stability and 
dynamics of small amplitude excitations of BEC arrays 
with positive scattering length against a uniform back- 
ground and in presence of a linear potential modelling 
the gravitational field. The linear stability analysis of 
the background for BEC arrays is investigated, and the 
region in parameter space where modulational instabil- 
ities develop is provided. We show that for particular 
parameters, shock wave propagation can also develop. In 
section V the main characteristic of Bloch oscillations for 
arrays with negative scattering length, are discussed. Fi- 
nally, in the last section VI the main results of the paper 
are summarized. 



As is well known, in the mean field approximation the 
wave function of a BEC in a trap potential V(r) satisfies 
the time-dependent Gross-Pitaevskii equation (CPE) 
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9*(r, t) 

at 



[-£-A + V(r) 
2m 



g \^(r,t)\ 2 ]^(r,t), (1) 



where go = 4Trh 2 a s /m, m is the atomic mass, a s is the 
s-wave scattering length which can be either positive (the 
case of 87 Rb atoms, with repulsive interactions, and with 
a s = 5.5 nm) or negative (the case of Li atoms, with 
attractive interactions and a s = —1.45 nm). The trap 
potential is assumed to be periodic along the z-direction 



V{r) = V{r + Xk), 



(2) 



where A is the spatial period of the potential and k is a 
unitary vector in the z-direction. A typical model poten- 
tial is 



V(r) — mgz + Vq(x, y) siu 2 (kz), 



(3) 



with experimental parameters A = 27r/fc ~ 850 nm, Vq ~ 
2.1h 2 k 2 /2m, and with an atomic population in each wells 
of N w 10 3 atoms Q. 

To make analytical studies it is convenient to look for 
solutions of Eq. (|l|) of the form 



AT/2 



n=-M/2 



(4) 



where the summation is performed over the number, Af, 
of minima of the trap potential which is assumed to be 
even. The functions $ ra (r) = $(r — r„) are assumed to 
be strongly localized around the site n of the potential 
and normalized to the mean number of atoms in the n-th 
well A = N/M w , 



$„(r)$„(r)dr = N 



(5) 



where N is the total number of atoms and Af w = Af + 1 
the total number of potential wells (hereafter the overbar 
will denote complex conjugation). This implies that the 
hopping integral 



\J, 



n,n+l | 



$ n (r)$ n+ i(r) dr 



-<€ N, 



(6) 



can be neglected, i.e. we assume J n , n +i ~ for all n. 
This is analogous to the well known tight binding ap- 
proximation of solid state physics |2lf |, valid for weakly 
overlapped condensates. Note that from (|5|) and (Q) the 
following normalization condition for the functions tp n (t) 
is obtained 
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or "finite density" 



Substituting (^) into ([!]) , multiplying by i>„ (r) and inte- 
grating over r, we readily obtain the system of coupled 
equations 

ihip n = E n ip n + U n \ip n \ 2 ip n - 

#n,f»-l^n-l - #n,n+l^n+l + 70™V>n, ( 8 ) 

where the overdot means time derivative, and 70 = m<?a 
with a = A/2 the distance between adjacent wells. In 
Eq. (@) we have denoted with 

= / £' V$ "I 2 + I^I^W] dr ( 9 ) 
the zero-point energy for the well n, with 

[4 = f||$„| 4 dr (10) 

the mean-field self-interaction energy, and with 

-2 



K, 



n,n±l 



[ r (Vf n V$ B± i) + (*n^*n±l)] * 

2m 



(11) 

the sum of the kinetic energy and the off-diagonal matrix 
element of the trap potentials between sites n and n + 1 
(i.e. the coupling constant between neighboring BECs). 
In the following we shall consider the case of equal con- 
stants E n = E, U n — U, K nn ±i — K. Introducing the 
new variables 

t-ft, ^ - (M) 1 ^*- 9 '*. (12) 

we can rewrite Eq. (||) in dimensionless form as 

+ ip n+ i + V'n-i - 2^n + 2o-|t/i„| 2 V'„ + 7 n- 0n = 0, 

(13) 



with a = — sign(a s ) and 7 = jq/K. Equation ( pL3| ) coin- 
cides with the well known discrete nonlinear Schrodinger 
equation |^], extensively investigated in the theory of 
nonlinear lattices. In the present context, however, be- 
cause of the normalization condition (^), the boundary 
conditions for Eq. (|l3|) cannot be choosen arbitrarily. In 
this paper we consider either periodic boundary condi- 
tion 



%j}(n) = ijj(n+j\f +1) 
in the case of finite lattices, or zero 
lim ip( n ) — 

1 71 — > OO 



(14) 
(15) 



lim ip(n) 

',— »±oo 



pe 



±1? 



(16) 



boundary conditions, in the case of infinite lattices (here 
p and $ are arbitrary constants). Eq. (|l3|) will then have 
also another integral of motion of the form 



(17) 



The mean field tight-binding model in Eq. (13) will be 
used in the following to study the dynamical properties 
of both ILM and ES excitations in ID BEC arrays. 



III. EXISTENCE OF INTRINSIC LOCALIZED 
MODES IN BEC ARRAYS 

As is well known, in the case of continuous ID BECs, 
soliton ground state solutions can exist. These were theo- 
retically predicted in p7| and experimentally observed in 
the case of positive scattering length (dark solitons) p8[ . 
Bright solitons can exist in BECs only when the effects 
of focusing nonlinearity balance the ones of effective dis- 
persion, this being possible for attractive inter-particle 
interactions or negative scattering lengths, only. Since 
bright solitons are very important for developing BEC 
applications (29), it is of interest to investigate their ex- 
istence also in BEC arrays. In this case, as we will see 
in the following, new possibilities can arise. In Ref. || 
bright solitons of the ES type were numerically investi- 
gated in ID BEC arrays with positive scattering lengths. 
Here we shall consider the case of bright solitons of ILM 
type in BEC arrays with both positive and negative scat- 
tering lengths. We fix 7 = , a = 1, in Eq. (|13|), i.e. we 
consider BECs in horizontal traps with negative scatter- 
ing lengths. ILM are then expected to exist when the 
coupling constant (which is responsible in the linear case 
of the spreading of the BEC wavefunction) is small in 
comparison with the nonlinear self-trapping interaction. 
From the normalization condition (|), and from defini- 
tions @, (|l|), one can estimate \U\/K = O(N), i.e. 
a stronger nonlinearity corresponds to a larger number 
of particles. For the sake of analytical developments is 
convenient to introduce the small parameter 




0(N-^ 2 ) < 1 



(18) 



representing the ratio between tunnel coupling and non- 
linearity, and recast Eq. (fT3J) it in the form: 



itp n + K(ip n+ t + ip n -x) + \<p n \ <Pn = 0, (19) 

with (p n = T/Kip n e 2l,t (the overdot here denotes the 
derivative with respect to r = t/n). 

In the following we discus ILM solutions of Eq. ( JT9| ) 
with different symmetry properties. 
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A. Symmetric ILM centered on a site 

Symmetric ILM solutions centered on a site (say n = 0) 
<P-n(r) = <p n {r), Mt)| = 1, (20) 

can be searched in the form 

oo oo 
ip n = G MT ^2 tfinjK 3 , U = ^2ujjK\ (21) 

j=n j=0 

with the frequency of the local oscillation oj equal for 
all sites. Substituting this expansion into Eq. (|l^) and 
performing straightforward algebra, we obtain 

uo = 1 + 2k 2 - 2k 6 + 2k 8 + 0(k 10 ), (22) 

with the site amplitude ip n satisfying the lattice equations 

<pi = [«- k 5 + n 7 + 0(n 9 )]e luJT , 
ip 2 = [k 2 - k 4 + k, 8 + 0{K W )\e iuT , 
^3 = [k 3 - 2n b + k 7 + 0{n 9 )]e luJ \ 
<p 4 = [k a - 3k 6 + 0{n 8 )]e luJT , 
<p 5 = [k 5 + 0{K 7 )]e^ T , 



from which we see that \(p n \ 3> l^n+il provided k«1. 

It is remarkable that the decay of the amplitude ip n is 
exponential as one move away from the n = site, and 
in the leading order can be approximated by 

tp n « exp(iwT - 77 |n|), (23) 

with r\ given by exp(— rj) = K << 1. Note that in accor- 
dance with our scaling, the limit n — > can be viewed 
as equivalent to the strongly nonlinear limit N — > 00 of 
the original physical system. Moreover, from Eq. ( |23| ) it 
follows that a bright onsite centered ILM can have for 
k = 0.1 about 91 % of all the atoms concentrated on its 
central site. For negative scattering lengths this implies 
that ILM can be stable only if N < N cr , where N cr is 
the critical threshold at which collapse phenomena occur 
(the effective coupling constant cannot be made smaller 
than a critical value, this being a feature of BEC arrays 
with respect to other nonlinear lattices). 

To check these results we have performed numerical in- 
tegrations of Eq. ([ll]) using as initial condition Eq. J23| ) 
with r = 0. We used a DOPRI8 integration routine 
p0| , based on a Runge-Kutta scheme of 7th-8th order 
with automatic stepsize control, so to combine speed with 
high precision (the relative errors was from 10~ 7 up to 
10 -13 ). The integration domain was taken large enough 
(800 sites) to avoid the influence of boundary conditions. 
The results are presented in Figs.[|a,b for the case, re- 
spectively, of weak and strong coupling among neigh- 
boring BECs in the array. We see that in the case of 



weak coupling the initial atomic population remains sta- 
ble for arbitrary long times (Fig. la), while for strong 
couplings, the initial distribution of atoms spreads out 
over the whole system (Fig. lb). 

FIG. 1. Distribution of the normalized atomic population 
according to numerical solution of Eq. ([L9|) for the case of 
weak coupling, k — 0.1, (a) and strong coupling, k = 0.8 
(b): this distribution corresponds to t — 10. The dashed line 
shows the ILM envelope given by Eq. (|||) . 

B. Anti-symmetric ILM centered on a site 

ILM of anti-symmetric type 

<A) = 0, ipi = 1, ip-„ = -ip n , (24) 

can be searched in the form (n > 1) 

00 00 
ipn = eliA)=T X! Pnjtf^ 1 , W = ^ UjK j . (25) 
j=n j=0 

Direct substitution into Eq. ( |l9|) gives 

uj = 1 + k 2 + k 4 + 2k 6 + 6k 8 + 0(k 9 ), (26) 

(p 2 - [k + k 3 + 2k 5 + 6k 7 + 0(K 9 )}e luJT , 
tp 3 = [n 2 + k a + 2k 6 + 5k s + O(n 10 )]e luJT , 
ifi = [k 3 + k 5 + k 7 + 0(n 9 )}e iuJT , 
<p 5 - [k 4 + k s + 0( K 8 )]e^ T , 



In analogy with the symmetric case, one can ensure 
that to leading order the shape of the ILM is well de- 
scribed by the function 

ip n w exp(iwr — 77 \n — 1|), n > 1. (27) 

Direct numerical simulations of Eq. ( p3|) with initial con- 
ditions given by Eq. (]27|)), showed a behavior of the 
atomic population distribution as a function of the cou- 
pling constant, similar to the one reported in Fig. [I]. 

Other types of bright symmetric and anti-symmetric 
ILMs, such as the ones which are centered between two 
sites, can also be constructed. A more detailed analysis 
of these excitations shows that the stability of ILMs de- 
pends on their symme try prope rty, as well as, on their 
localization extension |17]Jl9| , |3l| . Thus, symmetric on- 
site centered ILMs are stable but symmetric centered be- 
tween sites ILMs are unstable, while antisymmetric cen- 
tered between sites ILM are stable if k < n cr Rj 0.12, and 
unstable if n > n cr . 

For positive scattering lengths, different families of 
ILMs (as well as discrete fronts) representing dark ILM 
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(i.e. a constant atomic population along the array except 
for few sites in which the population is decreased), can 
also be found 0,11,11. 

It is important to note that, due to the invariance 
of Eq. ( p~9| ) under the transformation ip n — > (— l) n cp n , 
K — > — k, we have that the above bright (dark) ILMs, 
multiplied by a factor (—1)", are also solutions of Eq. jl9| ) 
for a s > ( a s < 0). This shows Eq. ( |Tj| ) possesses both 
bright and dark ILM solutions, for any sign of a s . 



IV. BEC ARRAYS WITH POSITIVE 
SCATTERING LENGTHS 



iq n + e 11 q n+ i + e %1 q n _i - 2q n + 2(p - \q n \ )q n = 0, 

(35) 

with boundary conditions 

q n — > pexp[ix(t)], at n — > ±oo. (36) 

In the next subsection we use this equation to study the 
stability of the background of a BEC array in presence 
of linear potentials (7 ^ 0). 



A. Linear stability analysis of the background 



From the experiment in Ref. |Tl| on vertical BEC ar- 
rays and from numerical simulations, it is known that 
the phenomenon of Bloch oscillations of ES can occur. 
In this section we shall use methods of nonlinear lattices 
theory j^,|25|,|34|,|35| to investigate Bloch oscillations of 
ES in BEC arrays with positive scattering lengths. For 
bright ES on zero backgrounds this was done in Ref. || , 
so we concentrate here only on the case of ES on top 
of non zero backgrounds. To this end it is convenient 
to introduce the background amplitude p directly into 
Eq. ( p"3| ) with a = — 1, this giving 

iip n + ip n+ i + <^„_i - 2ip n + (28) 
2{p 2 - |v?„| 2 )v?„ + jntfn = 

where cp = ip n ex^{—2ip 2 t). Let us consider the case of 
infinite number of lattice sites for which the integrals of 
motion in Eq. (pq) are written as 



and 



H = J2((^n-l~-p 2 )-^n\ 2 -p 2 f) 



(29) 



(30) 



Moreover, since the problem has nonzero boundary con- 
ditions we must have, for consistency with Eq. (p8|) [p5| , 
that 



— * pexp[i<I> n (t)] at n — ► ±00, 
with the phase $„(i) defined as 



$„(t) = n*y + *(*), 
' sin "ft 



X(t) = 2t 



It 



The following gauge transformation |3Jj3E 

a — id p~ l ~< nt 
allows to rewrite Eq. (p8|) as 



(31) 

(32) 
(33) 

(34) 



Let us consider a solution of Eq. ( |35| ) of the form 

q n = [p + a„(t)] exp[ix(t)], (37) 

where |a n (t)| <C p. Linearizing Eq. ( ^5|) around 
pexp[ix(t)] and expanding in Fourier modes 

00 

/3(M)= E a n(ty kn , 



n— — 00 



2tt 



«n(*) = 2^ / P(k 1 t)e- tkn dk 



we obtain 



(ft 



\k,t) = T k (t)\k,t), 



(38) 



where T& (t) is a 4 x 4 matrix whose nonzero elements are 



T, 



-,(12) _ _ T (21) _ 



(34) _ ^(43) _ 



k — M k 

Ti 



2[cos(7i — k) — cos(7i) — p 2 
2[cos(7i + k) - cos(jt) - p 2 



T (14) = T (23) _ ^(32) = T (41) = 2 9 
/c /c \z ' ' 



and with \k,t) a Bloch-Floquet state written in vector 
form 



\k,t) = 



( Pi(k,t) 
(3 2 (k,t) 

\ /%(-*, t) 



(39) 



(here /3(fc,f) = Pi{k,t) + i/3 2 (M) with 0j(h,t) real). 



FIG. 2. Floquet stability analysis of a k-modulated back- 
ground as a function of 7. Regions of stability are designated 
with "+". 
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Since is periodic we can use Floquet theory |3(| to 
study the stability. To this end we introduce the matrix 
Sfc(i) as a solution of the equation 



[p + p 2 a]e 



i[x+f"f>] 



dS k /dt = T k S k 



(40) 



satisfying the initial condition Sfe(0) = I (here I is the 
4x4 unit matrix) . It then follows that the background is 
unstable if |A,,(&:)| > 1, where \j(k) is the eigenvalue of 
Sfe(27r/7) (note that S k (t) is unimodular, detSfc(i) = 1, 
so that A1A2A3A4 = 1). In Fig. ^ the results of Flo- 
quet analysis, as obtained from numerical integrations of 
Eq. (f4(i|), are reported (in the numerical scheme we used 
a DOPRI8 procedure to integrate Eq. ( |4"c| ) and computed 
the spectrum of the matrix S& by reducing it to Hessen- 
berg form |}7) ) . The stability of the background was also 
investigated by direct numerical integrations of Eq. (|35| ) 
taking as initial conditions a uniform background of am- 
plitude p — 1 modulated by a sine wave of wavenumbcr 
k n = 2im/L (we used L = 400 lattice sites). 

FIG. 3. Modulational stability of the background at 7 = 0, 
for an initial modulation with k ~ 0.31 



FIG. 4. Modulational instability at 7 = 0.5, for an initial 
modulation with k ~ 0.31 

These results are reported in Figs. ||, |] for the cases 
7 = 0, and 7 = 0.5, respectively. We see that in pres- 
ence of the linear potential a modulational instability can 
develop, in agreement with the Floquet theory analysis 
of Fig. |^. A good agreement was found also for other 
values of 7 and other initial conditions (at 7 = the sta- 
bility was checked only by numerical integrations of the 
original equation, since in this case t = 2tt/j — > 00 and 
the Floquet theory becomes unapplicable) . In particu- 
lar we checked that the modulational instability of the 
background develops at later values of t, as 7 increases, 
and that an initial background with k = tt becomes again 
stable for 7 ~ 5.4, in agreement with Fig.@. 



B. Dynamics of small amplitude pulses 



The fact that Eq. fl28p possesses an integral of mo- 
tion of type (p9|), implies that it cannot support solitary 
waves |23|1 , i.e. it can not have solutions moving with 
constant velocity v and depending on n, t, through the 
combination n — vt (this is true also for 7 = 0). On the 
other hand, one can show that localized smooth excita- 
tions of small amplitude can propagate along the array 
with a relatively weak distortion. To investigate this dy- 
namics, it is convenient to introduce a small parameter 
p -C 1 defined as the square root of the ratio between 
the deviation from the background and the background 
itself, and consider the case of very small 7, i.e. we take 
7 = o(p 3 ). This allow us to look for solutions of the form 



(41) 
0(p 4 ) 



where a = a,Q + pfa\ + 0(p ) and <f> — cf>o + p <f>\ 
are real functions of the slow variables £ = pn, T = pt, 
t = p 3 t, considered to be continuous and independent. 
We can then perform a multiple scale expansion |2l| of 
Eq. (j35|), substituting jt with St, with 6 ~ o(l). A 
straightforward algebra shows that the expansion equa- 
tion at the hrst order in p is identically satisfied by the 
substitution (|4l|). 

FIG. 5. Time dependence of coefficients of Eq. (^) for 
5 = 0.1, p=l. 

The equations at the orders 0(p 2 ) and 0(p 3 ) are 
1 



a 



4p 



[d T 4>o + 2sin(<5r)9 5 0o], 



(42) 



and 

9t4>o + 4sin(<5T)cV<9£</>o + 

4[(sin(,5T)) 2 - p 2 cos(5r)]<9f 0o = 0, (43) 

respectively. It is worth to note that Eq. (^) at times 
5tq = ±7r/2 changes from hyperbolic (supporting wave 
propagation if cos((5r) > 0) to elliptic. Let us consider 
the dynamics of the small amplitude excitation during 
the time in which Eq. (^) is of the hyperbolic type. We 
introduce a new running variable £ — £ — V(t)T where 



V(t) = 2sin((5r) + 2 P v / cos{6t) 



(44) 



can be interpreted as slowly varying velocity of the wave 
packet (for the sake of definiteness we have chosen only 
one branch of the solution, the other branch being char- 
acterized by the "velocity" with opposite sign). One can 
then show that Eq. ( |43| ) is satisfied for arbitrary pair of 
functions oo((,t), <po(Ci T ) linked by 



1 



ao{(i T ) = -\/cos((5t)<9 C (/>o(C,t). 



(45) 



In order to find the dependence of those functions on r, 
one has to consider the equations of the forth and fifth 
orders in p (more precisely the condition of their compat- 
ibility). After tedious but straightforward calculations, 
one arrives at the following Korteweg - de Vries (KdV) 
equation 

d T a + 6i(r)a 9 c a + b 2 {T)d 3 a + fro( r ) a o = 0, (46) 

with slowly (if S = o(l)) varying coefficients bj(r) given 
by 

, , , 1 , , , , 3[cos((5t)] 3 / 2 -psmtSr) 

b r = 7^ tan {St , b x r = 2^ V U , F —, 

4 cos(or) 

h(r) = — (4psm(ST) + p 2 \cos(St)}^ 2 - 3[cos(<5r)] 3 / ; 
12p V 

FIG. 6. Evolution of the localized excitation in a BEC ar- 
ray, governed by Eq. (fl6|). The initial wave profile corre- 
sponds to the pure KdV soliton (Eq. @) with parameters 
UQ = 1, (3 = 1. 
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In Fig.|6| we report the time evolution of a numerical 
solution of Eq. (|4q ) with initial condition 

a (C,T) = 3 V sech 2 (^^), l c = U^j 1 . 

Two important features can be seen from this figure. The 
first is the existence of a breaking time 7V below which 
the profile of the wave is distorted but its evolution is 
smooth (in Fig. ^ Tbr ~ 0.5 ). At t = tv the breaking of 
the wavefront occurs, after which a train of solitary pulses 
(each of them being a KdV soliton) is emitted. This is 
similar to the generation of shock waves in a fluid (non- 
linear Schrodinger lattice shock waves were first studied 
in Ref. plj). Note that at the breaking time the coef- 
ficient b^jbr) is exactly zero (see Fig. ||), and Eq. J46| ) 
reduces to a known shock dynamics equation p8| ] . An es- 
sential difference between the case at hand and the case 
7 = treated in Ref. pi] ], however, is that in this last 
shock waves can appear only for a given carrier wave 
background, while here the linear potential (gravitational 
field of BEC arrays) makes them possible to exist for 
any wavenumber. Note that for the observation of shock 
waves it is important that 8 is small compared with the 
amplitude of the background, this ensuring the time at 
which the group velocity dispersion becomes negligible, 
be long enough for shock waves to develop. 

The second important feature which emerges from 
Fig. ^| is the possibility to have bright pulses propagating 
on top of a nonzero background (i.e. on top of a constant 
atomic population along the array). Thesepulses actu- 
ally are ES, similar to those predicted in [g0| for 7 = 0. 
For 7^0, these excitations can undergo Bloch oscilla- 
tions. To show this we have performed numerical inte- 
grations of Eq. (|3|), with a sech initial profile of large 
amplitude (i.e. comparable with the background level). 
The results are reported in Fig.J^ (the numerical scheme 
used is the same of the previous section). We note that 
for the same parameter values, the bright excitation (i.e. 
above background) breaks down while the dark one (be- 
low background) remains stable for long times. 

FIG. 7. Bloch oscillations of bright (a) and dark (b) local- 
ized excitations. Initial conditions: 

9n(0) = p ± 2f?Sech[2f7n]exp[#„(0)], p = 0.1, 77 = 0.04, 
7 = 0.1, 0„(O) = -fn. 

An important conclusion following from these numeri- 
cal studies is that bright and dark solitons of KdV type 
on top of homogeneous backgrounds, can exist in BEC 
array with positive scattering lengths. 

V. ARRAY OF BEC WITH NEGATIVE 
SCATTERING LENGTHS 

A distinctive feature of BECs with attractive inter- 
actions is the possibility of collapse when the number 



of atoms in the condensate exceeds a critical value N c . 
From GPE one can predict N c ~ 1400 for 7 Li atoms, a 
value which was confirmed experimentally in Ref. |39|. 
Arrays of optical traps can be used for manipulation of 
BECs with negative scattering length, if their wells are 
loaded with number of atoms less than this critical value. 
In this situation it is reasonable to consider the dynam- 
ics of ES on zero background, as described by Eq. ( |l~3| ) 
with <t = 1. Bloch oscillations of bright ES were numeri- 
cally investigated in Ref. |m| . From these simulations, a 
strong dependence of the dynamics on the amplitude of 
the initial wave was found (large amplitude wavepackets 
are quickly reduced to fragments, while small amplitude 
excitations keep their integrity over many oscillation pe- 
riods). From numerical studies is also known that the 
amplitude of a dynamically localized ES can oscillate in 
time. Using the method of the previous section, we can 
develop analytical considerations on Bloch oscillations in 
BEC arrays with negative scattering lengths which ex- 
plain the origin of these oscillations. We consider the 
case of small amplitudes and small 7, i.e. we assume, as 
before, /1 = y/j to be a small parameter of the problem. 
Introducing the variable 

q n = i>n exp(-i7nt), (47) 

Eq. (|l~3| ) can be written as 

iq n + cos(jt)(q n+1 + (?„_i - 2q n ) + 
isin( 7 t)(g„ + i - q n _!) + 2\q n \ 2 q n = 0, (48) 

whose solution can be searched of the form 

<7n = M<2(£,r), (49) 
with £(t) = /in — x(t), t = /i 2 t, and 
2 

x = -(1-cost). (50) 

Note that x(t) = 0(/i 3 ) = 0(7 3/2 ) and thus the substi- 
tution is self-consistent (at /i 3 i <C 1 the term cos(r) can 
be expanded in Taylor series). After calculations ana- 
logues to those of the previous section, one arrives at 
the following nonlinear Schrodinger (NLS) equation with 
periodically varying dispersion 

id T Q + cos(T)dlQ + 2\Q\ 2 Q = 0. (51) 

The dynamics of this equation was numerically investi- 
gated starting with initial conditions of the form 

Q(£,0)=2»pech(2»;(£-G,)). (52) 

For the numerical code we used a split-step fast Fourier 
transform technique [fl0| with a time step At = 0.001 
and a space step A£ = 0.02 (in normalized units), corre- 
sponding to 1024 grid points in the discretization domain 
taken as -10 -r 10 (absorbing boundary conditions were 
used to simulate an infinite domain). The accuracy of 
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the scheme was checked by monitoring the conservation 
of the number of atoms and of the Hamiltonian, which 
was within ± 0.1 % for all runs. The results are depicted 
in Fig.| 

FIG. 8. (a) Evolution of the localized excitation governed 
by Eq. j5l|), for 77 = 0.5. and (b) Decay of its amplitude. 

From Fig. ||a) we see that while the ES is executing 
spatial oscillations, its amplitude is oscillating in time. 
In Fig. ||b) we report the time dependence of the cen- 
ter of the wave on a longer time scale, from which we 
see that the amplitude while oscillating is also decaying 
(note, however, that this dynamics reproduces the dy- 
namics of a real array only for times r such that At < 1). 
A more detailed numerical investigation shows that the 
amplitude can be either decaying or growing in time, de- 
pending on initial conditions (the growth can reach also 
the 60% of the initial amplitude of the pulse) . This im- 
plies that if at i — the number of atoms in one well is 
below the critical value No(t = 0) < N cr for collapse, the 
growth in amplitude due to the linear potential, can in- 
duce collapse after a time t — t cr at which No(t cr ) = N c . 

VI. CONCLUSIONS 

The dynamical properties of BEC arrays with positive 
and negative scattering lengths have been studied in the 
framework of a discrete nonlinear Schrodinger equation 
derived from the mean field GPE with a tight-binding 
approximation. The interplay between macroscopic in- 
ter site tunnelling and nonlinear self-trapping was shown 
to be responsible for the appearance of different types 
of DMS. In particular we showed, both analytically and 
numerically, that for a small ratio between the tunnel 
coupling constant and the nonlinear interatomic inter- 
actions, ILM solutions can exist. For BEC arrays with 
non-zero backgrounds the modulational stability problem 
was investigated and the existence of bright and dark ES 
was discussed. The problem of Bloch oscillations of en- 
velope solitons in arrays with positive scattering lengths 
was also analytically and numerically investigated. We 
showed that at the lower orders of a perturbative expan- 
sion, the dynamics of small amplitude excitation evolves 
according to a KdV equation with time dependent coef- 
ficients. In this case the possibility of shock waves for- 
mation, in BEC arrays with positive scattering lengths, 
was explicitly displayed. In the case of arrays with neg- 
ative scattering lengths, the dynamics of small ampli- 
tude excitations was described in terms of a nonlinear 
Schrodinger equation with a periodically varying disper- 
sion. This equation was used to show the presence of 
amplitude oscillations during Bloch oscillations, as well 
as decay or growth of the excitations, depending on the 
initial conditions. The results of this paper clearly show 
the complexity and the wide range of behaviors which 
can arise in the system. 



In closing this paper we feel compelled to discuss to 
which extend the phenomena presented in this paper 
could be observable in real BEC arrays. Our model is 
based on a tight binding approximation, this putting re- 
strictions on the shape of the wavefunctions, as well as, 
on the potential profile, to be used. Although a detailed 
analysis of the experimental settings for which the tight 
binding approximation is accurate has not yet been done, 
there are experimental situations in which this is obvi- 
ously true. Thus, for example, if the potential barrier is 
wide and tall enough to reduce the tunnelling probability 
among adjacent BECs, the overlapping of the wavefunc- 
tions is certainly small and the model is accurate. In this 
case DMS of the type discussed above should appear. 

In particular, the possibility to observe ILMs in real 
experiments should not be overlooked. We remark that 
from Eqs. (0), (|l2|), it follows that the peak amplitude, 
i.e. max(cj)n(t)) , is of order 0(1), so that the atomic 
population can be localized on very few sites of the array. 
In the experiment in Ref . [ pd[ , the peak density was hq = 
10 13 /cm 3 , the mean field energy gono » ks • 4nK (fee 
being the Boltzmann constant), this giving the parameter 
k in Eq. ( O ) to be k w 1. To observe ILM one needs to 
have at least k ~ 0.1. This can be achieved either by 
increasing the number of atoms in each well (from 10 3 
to 10 4 ), or by reducing the tunnelling constant K, which 
exponentially depends on the lattice constant (distance 
between potential wells). When the number of atoms 
in the wells is increased, however, a loss of coherence of 
the condensate can occur, a phenomenon observed in the 
experiments reported in Ref. |3^ ], On the other hand, 
k can be reduced also by changing the atomic scattering 
length a s using the Feshbach resonances p3| ], so that it 
should be possible to find experimental settings for which 
ILM can be observed. 

As to ES on finite backgrounds discussed in Sec- 
tion IV, we remark that the number of atoms in- 
volved in these KdV solitons can be estimated as N s = 
^HP X^o a (C/^c)^C = For the parameters used 

in the numerical simulations of Fig. ^, with S = /x 3 = 
0.1, Mo = /3 = 1 we have fj, w 0.46, l c = 3.464, and for an 
array of Rb atoms with 1000 atoms in each well, we have 
N s ~ 6370 atoms. These bright (dark) solitons on top 
of a background, could be created by using a laser beam 
applied for short time to a uniform BEC array, so to cre- 
ate a local enhancement (depletion) of the potential. In 
this case the dependence of the soliton velocity on its am- 
plitude is an interesting parameter to measure since, in 
contrast with nonlinear Schrodinger solitons, the velocity 
of KdV solitons depends on the amplitude of the wave 
(number of the atoms in the condensate). This could 
provide a way to experimentally check our results. 

We hope that experiments in this direction will be soon 
performed. 
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